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We have used the Path Integral Monte Carlo method to simulate a monolayer of molecular hy- 
drogen on graphite above 1/3 submonolayer coverage. We find that at low temperature and as the 
coverage increases the system undergoes a series of transformations starting from the \/3 x \/3 com- 
mensurate solid near 1/3 coverage. First, a phase is formed which is characterized by a uniaxially 
compressed incommensurate solid with additional mass density modulations along the same direc- 
tion which can be viewed as an ordered domain-wall solid with a characteristic domain-wall distance 
which depends on the surface coverage. At low temperature and higher coverage there is a transition 
to an incommensurate solid which is rotated relative to the substrate commensurate lattice. As a 
function of temperature the domain-wall ordering first melts into a fluid of domain-walls and at 
higher temperature the solid melts into a uniform fluid. Regardless of the large amplitude of quan- 
tum fluctuations, these phase transitions are analogous to those in classical monolayer fllms. Our 
calculated values for the surface coverage and temperature where these transitions occur, the cal- 
culated structure factors and specific heat are in general agreement with the available experimental 
results with no adjustable parameters. 

PACS numbers: 64.60.Fr, 67.40.-w, 67.40.Kh 



I. INTRODUCTION 

The commensurate-incommensurate transition has 
been expensively studied in classical monolayer films 
both theoretically and experimentally. Monolayers of hy- 
drogen or helium are quantum mechanical systems and, 
in principle, one might suspect a different behavior due 
to the effect of strong quantum fluctuations. 

The monolayer phase diagram of molecular para- 
hydrogen adsorbed on graphite as inferred from experi- 
mental studies 01 is shown in Fig. |l| This phase diagram 
was originally drawn based on the anomalies found in the 
specific heat|^, ^, ^; for the characterization of the vari- 
ous phases, low energy electron diffraction (LEED) |^ 
as well as neutron diffraction ^ ^ studies have been 
carried out. At 1/3 coverage the molecules condense on 
the surface of graphite in a -v/S x a/3 commensurate solid. 
In the low coverage region (p < 0.6) (in units of the 
y/3 X a/S commensurate density), it forms a commensu- 
rate solid-gas coexistence phase at low temperature and a 
2D gas phase at higher temperature (T > lOK). For cov- 
erages 0.6^/5^0.9, as a function of temperature, there is a 
transition from a commensurate solid cluster phase to a 
commensurate solid phase with vacancies at higher tem- 
perature and to a 2D gas at even higher temperature. At 
density somewhat higher than unity the so-called a-phase 
is formed which is believed to be a striped domain-wall 
solid phase and at higher temperature it transforms to 
the so-called /?-phase and to a fluid phase at even higher 
temperature. At high densities a transition take place to 
a triangular incommensurate solid phase. 



The phase diagram of molecules formed from the iso- 
topes of the hydrogen such D2 and HD molecules phys- 
iorbed on graphite has also been studied ||l^, |ll| and it is 
similar to that of H2 on graphite but more complex. 

We have recently studied the phase diagram of 
molecular hydrogen on graphite at and below 1/3 cov- 
erage using Path Integral Monte Carlo (PIMC) simula- 
tion. Our computed phase diagram was in general agree- 
ment with that inferred from the experimental studies 
for this coverage range. In this paper we extend the 
PIMC calculation to study the phase diagram above the 
commensurate density which is less well understood both 
theoretically and experimentally. At coverages above 
p ~ 1.05 the monolayer undergoes a commensurate- 
incommensurate(C-IC) transition. While the existence 
of the so-called a and /3 phases were known from the 
specific heat anomalies, their characterization came from 
LEEDH, and neutron scattering! l) studies. Using 
LEED experiments Cui and Fain|^, |7 observed a uniax- 
ial IC solid with striped superheavy domain walls and a 
rotated triangular IC solid at higher coverages p > 1.23. 
Freimuth, Wiechert, and Lauter||lj (FWL) presented a 
neutron diffraction study of the C-IC transition and their 
results are in agreement with the LEED results. They 
examined the commensurate-incommensurate solid tran- 
sition, especially the striped domain-wall solid phase. In 
the Qf-phase the diffraction intensity has a main peak 
characteristic of a compressed lattice and a satellite peak 
on the lower side of the commensurate peak wave- vector 
position kc — 1.703A~^ that arise from the spatial mod- 
ulation due to ordered striped-domain walls. As cover- 
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age increases, the separation of the two peaks increases 
and the sateUite intensity decreases. At higher tem- 
peratures the peak height drops and the satelhte peak 
vanishes which impUes that the commensurate domains 
vanish and the system becomes an isotropic fluid phase. 
As coverage increases, the first molecular hydrogen layer 
forms a rotated incommensurate solid (RIC) phase. Neu- 
tron diffraction studies show a sharp and intense peak at 
wave-vector k = 1.97A~^ at the density p = f.34 which 
reveals an IC equilateral triangular structure. This phase 
is continuously compressed as the coverage increases up 
to the highest density allowed before layer promotion. 
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FIG. 1: Phase diagram of molecular hydrogen adsorbed 
on graphite from Ref. reproduced here for easy refer- 
ence. Density 1.0 corresponds to the complete commensurate 
\/3 X \/3 phase. 

In parallel to these experimental investigations several 
important theoretical studies were undertaken to under- 
stand the commensurate-incommensurate transition! 13 
[T^ , [isl , [l^ , in physiorbed systems. The phases and 
the phase transitions which were theoretically predicted 
and studied can be understood as classical phenomena 
due to the strain induced on the adsorbate molecules of 
the monolayer film due to the substrate periodic struc- 
ture which arise from the competition between adsorbate- 
adsorbate and adsorbate-substrate interactions. Molecu- 
lar hydrogen and atomic helium physisorbed on graphite 
are quantum films characterized by strong zero point mo- 
tion. Therefore, one could question the degree to which 



a classical picture might be valid and might expect new 
phenomena and phases to occur. Motivated by these 
thoughts we extended our earlier investigation JT2[ to 
study this system above 1 /3 coverage up to layer comple- 
tion starting from the known hydrogen molecule-molecule 
and hygrogen-graphite interactions |l], ^ 23 and 
using the PIMCjl^ which is a Quantum Monte Carlo 
technique adopted for the study of strongly interacting 
quantum films by Pierce and Manousakis . 

We have simulated the first layer of molecular hydro- 
gen on graphite with a variety of simulation cells that 
are appropriate for examining different phase regions of 
the phase diagram. We have applied periodic boundary 
conditions along both directions parallel to the surface of 
the substrate. We have computed expectation values for 
the total energy, the static structure factor, the proba- 
bility distribution, and the specific heat by means of the 
path-integral Monte Carlo simulation method using the 
multi-level metropolis method. In the multi-level bisec- 
tion method we have used level 3 which optimizes the 
acceptance ratio. To thermalize the system, we typically 
carried out of the order of 15,000 MC steps and we car- 
ried out of the order of 2,000-20,000 MC steps in order 
to compute observables. 



II. UNIAXIALLY COMPRESSED SOLID WITH 
ORDERED DOMAIN WALLS 



In this section we discuss the region above the com- 
mensurate coverage where experimentally the ordered 
domain-wall solid was found. We carried out a simu- 
lation at two densities in this region, the same densities 
discussed by FWL, where they inferred the existence of 
relaxed superheavy domain walls. 

We used simulation cells that can accommodate the 
structures that have been proposed by FWL. a) First, 
we consider p — 0.0694A~^. We have chosen a simula- 
tion cell with dimensions x— 5\/3 agr and y= 22 agr, 
where agr = 2.459A is the carbon-carbon distance on the 
graphite surface and 80 hydrogen molecules to produce 
the above density. In Fig. || we give the contour plot 
of the calculated probability distribution where we see 
that two commensurate-solid domains are separated by 
the incommensurate solid domains. The solid structure 
is uniaxially compressed along our y direction such that a 
new row of molecules is introduced for every 8 rows. No- 
tice that the period in the cc-direction is -s/Sogr and there 
is modulation along the y direction with period llugr, so 
that the wavelength of this striped domain- wall mod- 
ulation is Xs = 27.049. b) Second, we have performed 
the simulation at the density p — 0.0716A~^. We have 
chosen a simulation cell with dimensions x = 5\/3agr and 
y — IGugr and 60 hydrogen molecules in order to achieve 
the above mentioned density. The calculated contour 
plot of the probability distribution is shown in Fig. |^. 
Notice that in this case also there is an ordered striped 
domain-wall solid phase along our y direction. Namely, 
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FIG. 3: Contour plot of the probability distribution at the 
density 0.07161"^ for r= 1.33 K. The simulation cell is 
21.295 A X 39.344 A (60 molecules). The filled circles indicate 
graphite adsorption sites. The solid structure is uniaxially 
compressed along our y direction. The striped domain-wall 
solid phase can be seen. 



FIG. 2: Contour plot of the probability distribution at the 
density 0.0694 A"^ and r= 1.33 K. The simulation cell is 
21.295 A X 54.098 A (80 H2 molecules). The filled circles 
indicate graphite adsorption sites. The solid structure is uni- 
axially compressed along our y direction. In addition, two 
commensurate solid domains separated by the incommensu- 
rate solid domains can be seen which implies a mass density 
wave along the same direction. 



the amplitude of the molecular density wave is modulated 
along the y axis with a period of about half our cell size 
along the y axis. Notice that there are two commensurate 
solid domains separated by denser regions. 

Our computed static structure factor S{k) for p — 
0.0694A-2 is shown in Fig. |. The main Bragg peaks of 
this structure occur at ki = (1.475 A~\0.929 A"^) and 
k2 = (0,1.858 A~^). The satellite peak occurs at ^'sat= 
(0,1.626 A^^). The experimental values of the magni- 
tude of the wave- vectors at the peaks are k'^ain— 1-743 
and k'^^f= 1.632 A~^, which compare well with our 
computed values of fcf= 1.743 A~^ and k'^^^= 1.626 
We believe that FWL could not observe the peak at k2 
because of the strong interference with the (002) graphite 
peak. 

The interpretation of these results is as follows: Ana- 



lyzing the contour plot of Fig. g we find that a) the solid 
appears uniaxially compressed along the y direction, by 
adding another row for every 10 rows of molecules and 
spreading them evenly, while along the x direction the av- 
erage spacing between the molecules remains that of the 
commensurate solid, b) superimposed there is a density 
modulation along the y direction which has wavelength 
several times greater than the average nearest neighbor 
distance but of small amplitude. The actual size of the 
wavelength As of this striped domain-wall modulation is 
half of the length of our simulation cell in the y-direction, 
i.e., Xs = 27.O5A, as discussed in the previous para- 
graph. The two main Bragg peaks correspond to the 
unit vectors which span the reciprocal space of this uni- 
axially compressed triangular solid structure. The satel- 
lite peaks are due to lateral modulation in our y direction 
and are located at hsat = ^1,2 — kg. ki_2 correspond to 
the wave vector of each of the two main Bragg peaks and 
kg = (0, 27r/As), where As is the wave length of a unit cell 
for the striped domain-wall solid structure in the modu- 
lated direction. Using this As value, k2 = (0, 1.858A^^) 
and ks = (0, 0.232A~^) (obtained by using the values of 
As — 27.O49A found by inspection of Fig. 0) we can find 
that the satellite peak position is at ksat = (0, 1.626A^^), 
which agrees well with our peak value. 

The structure factor for the case of p = 0.0716A~^ 
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FIG. 4: Static structure factor S{k) at the density 0.0694A ^ 
(and corresponds to Fig. ^) for T— 1.33 7^ as a function of 
the magnitude of k. The main Bragg peaks are at ki= 1.743 
and k2= 1.858 . There is a satelhte peak located at 
1.626 A-\ 



30 



CO 




1.5 1.6 1.7 1.8 1.9 2 

k(A-') 

FIG. 5; Static structure factor S{k) at the density 0.0716A~^ 
(and corresponds to Fig. ^ for T— 1.33 K. The main Bragg 
peaks are at ki= 1.759 and k2— 1.916 A~^. Notice that 
the sateUite peak at k'^at~ 1-597 A~^ is significantly weaker 
than the satelhte peak at the lower density shown in Fig. 0. 



(Fig. |) is shown in Fig. |. The main peaks are ki = 
(1.475 A-\0.958 A-i) and ka = (0,1.916 which 
correspond to the uniaxially compressed solid produced 
by adding another row of molecules for every 6 rows of 
molecules and spreading them evenly. This, in addition, 
forms a superscructure with wavelength As = 19.672A 
which gives a satellite peak at k^„(= (0,1.597 1-^). No- 
tice that the scattering intensity, relative to the previ- 
ously discussed case of p = 0.0694A~^, has the peaks 
shifted at higher values of k as expected and the intensity 
of the satellite peaks is decreased as in the experiments. 



III. DOMAIN WALL MELTING 



In Fig. g we show the contour plot of the the probabil- 
ity distribution where a striped domain-wall fluid phase 
(at T = ll.llA') and a fluid phase (at T = 16.67) are evi- 
dent. Notice that the stripe domain- wall fluid phase (left 
part of Fig. |^) is characterized by mobile commensurate 
and incommensurate domains. To understand this con- 
tour plot we need to be reminded of the periodic bound- 
ary conditions used in our simulation and the fact that 
these domains cannot intersect each other. This implies 
that one domain will oscillate back and forth between 
the neighboring domains. We can determine the melt- 
ing temperature of the domain- wall solid phase from the 
temperature dependence of the specific heat and from the 
temperature dependence of the peak height of the static 
structure factor S{k). 




FIG. 6: Contour plot of the probability distribution at the 
density 0.0716A~^ at two temperatures 11.11 K and 16.67 
K. 



The calculated specific heat as a function of tempera- 
ture is shown in Fig. |^ and is characterized by two peaks, 
the first corresponds to the melting of the domain-wall 
solid while the second peak indicates the melting of the 
stripes and the sohd into a fluid. Our computed values 
for the melting temperature are generally lower that the 
experimental values. For example, for p ~ O.G716A~^ 
we find that at T = 9.{)%K the striped domain- wall 
solid phase undergoes a transition to the domain-wall 
fluid phase and at T = 12.5K the /3-phase becomes an 
isotropic fluid phase. Factors for obtaining lower values 
for the critical temperature than the experimental val- 
ues could be the finite-size effects and the interaction 
strength used to describe the hydrogen-graphite interac- 
tion. We also computed the static structure factor shown 
in Fig. ^ and we studied the temperature dependence of 
the height of its first main peak which is shown in Fig. ||. 
The peak height significantly decreases near the melt- 
ing temperature of the domain-wall solid and near the 
domain-wall evaporation temperature. 
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versus temperature at the density 0.0716A 
cell is 21.295 A x 39.344 A (60 molecules). 



The simulation 



FIG. 7: Specific heat versus temperature at the density FIG. 9: The peak height of the static structure factor S{k) 
0.0716A~^. The long-dashed line is a spline fit to the spe- 
cific heat values and is a guide to the eye. The filled circles 
are the experimental specific heat at the density 0.0716j4~^. 
The simulation cell is 29.813 A x 39.344 1 (84 H2 molecules). 
Our computed values for the melting and evaporation temper- 
ature are the peak positions, 9.09 K and 12.5 K, respectively. 
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FIG. 8: Static structure factor S{k) at the density 0.0716A"^ 
for various temperatures. 



vaco and McTague[|13| . The angle of the incommensurate 
lattice relative to the graphite lattice is approximately 
5° in good agreement with the experimental valueQ for 
the above density. The calculated static structure factor 
(shown in Fig. |ll|) has one sharp peak at k= 1.99 cor- 
responding to an unregistered equilateral triangular lat- 
tice in agreement with the value of k= 1.97 reported 
from the neutron diffraction experimental study|l|. 




FIG. 10: Contour plot of the probability distribution at the 
density 0.08491"^ The simulation cell is 38.331 A x 22.131 
A {72 H2 molecules). 



IV. ROTATED INCOMMENSURATE SOLID 



The first layer of molecular hydrogen adsorbed on 
graphite forms an incommensurate solid phase before the 
first layer is complete. We have carried out a simulation 
at the density p = 0.0849A~^ using a simulation cell 
with dimensions x= 9 \/3 agr and y= 9 agr- The calcu- 
lated probability distribution(Fig. |l^) also clearly shows 
the presence of an equilateral triangular solid structure 
which is not registered with the underlying graphite lat- 
tice and it is rotated relative to the graphite lattice. 
This rotation was first predicted and discussed by No- 



V. ENERGY CALCULATION 

In Fig. |l2| we compare the energy per hydrogen 
molecule as a function of the surface coverage for the 
cases where the full potential is used (solid line) and 
the case where only the laterally averaged potential is 
used (dashed line). There are several important observa- 
tions which will make based on the results presented on 
the this figure. First the curve obtained using the lat- 
erally averaged potential has a minimum at the density 
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FIG. 11: Static structure factor S{k) at the density 
0.0849A"^. The simulation cell is 38.331 A x 22.131 A (72 
H2 molecules). The sharp peak position is k= 1.99 A~^. 



FIG. 12: The dashed line is the energy per hydrogen molecule 
on the graphite substrate using the laterally averaged poten- 
tial as a function of surface coverage. The solid line is the 
same quantity when corrugations are included. 



VI. CONCLUSIONS 



Pq = 0.0705A~^. At this equilibrium density the system 
forms a triangular lattice which fills the entire system and 
below this density the system undergoes a phase separa- 
tion and forms a solid-vapor coexistence phase. There- 
fore below Pq the energy as a function of density is higher 
for a finite-size system or for a system which is forced to 
be uniform at this lower density. Above this density the 
system is a compressed triangular solid until the promo- 
tion density|12 . The solid line corresponds to the en- 



ergy per molecule when we used the hydrogen-molecule- 
graphite interaction potential with the corrugations. The 
minimum in this case occurs at pc = 0.0636A~^ which 
corresponds to the VS x ^/S solid. Below that density the 
system is unstable to formation of solid clusters and this 
part of the phase diagram was investigated in Ref . . 
Above this density up to a density one can notice that 
the presence of the substrate corrugations moves the min- 
imum from po at much higher energy but there is a fea- 
ture at the same density produced by the H2-H2 inter- 
action. The value po of the equilibrium density for a 
smooth substrate plays a significant role even when the 
full potential with the substrate corrugations is used be- 
cause it is the density which minimizes the energy due 
to H2 molecule-molecule interaction. Clearly due to this 
term which has a preference for a higher density than the 
commensurate density, the solid line is not convex which 
a sign of instability of a uniform phase. These results are 
obtained on lattices which frustrate the proper periodic- 
ity of the domain walls. If the calculations arc done on 
lattice which can accommodate the periodicity of these 
superlattice structure the energy is lowered. As an exam- 
ple we also give (solid circle) the energy obtained for the 
stripe state of Fig. |3[ Notice that the energy falls below 
that obtained with choice of lattices unfavorable for this 
ordered domain-wall state. 



Our studies of H2 monolayer on graphite above 1/3 
coverage clarifies the nature of the phase diagram. The 
boundaries of the phase diagram were well-known from 
the specific anomalies. We find that the phase near 1/3 
coverage is the known V3 x commensurate solid. At 
higher densities the monolayer undergoes a transition to a 
uniaxially compressed solid with density wave stripe-like 
modulations along the same direction. This phase was 
first identified by LEED[|j studies; while the subsequent 
neutron diffraction studies confirm the existence of or- 
dered domains walls, the peak which is the signature of 
the uniaxial compression is not seen because of the (002) 
reflection from the underlying graphite lattice. In addi- 
tion in both studies the domain walls are assumed to be 
superheavy. As it is clear from both our contour plots 
and our structure factor the walls are not significantly 
denser than the rest of the system; we find that the addi- 
tional molecules added to the monolayer to increase the 
coverage above the commensurate coverage are used pri- 
marily to uniaxially compressed the solid and they are 
only partially used to construct the walls. Our calcu- 
lated specific heat for this coverage range as a function 
of temperature has two anomalies at two temperature 
values. From the calculated contour plots and the struc- 
ture factor, the lower temperature anomaly is identified 
as the transition from the ordered domain-wall solid to a 
solid where the domain walls form a fluid and the higher 
temperature anomaly is caused by the total melting of 
the solid to a uniform fluid. 

At even higher density we flnd that the uniaxially com- 
pressed solid with the ordered domain walls undergoes a 
transition to a triangular solid which is incommensurate 
with respect to the substrate lattice and in addition the 
triangular solid of the adlayer is rotated by an angle of 
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the order of 5° with respect to the underlying substrate 
lattice. This is in complete agreement with both the the- 
oretical prediction and experimental findings (l[ ||. 
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